library(tidyverse) #ggplot2, dplyr, tibble, etc.
library(plotly) #use ggplotly() for interactive plots with scoll-over IDs
library(knitr) #use kable() to make formatted tables
library(Rtsne)
library(glmmTMB)
library(lme4)
library(sjPlot)
library(ggpubr) #ggarrange() figure wraps
library(plyr)
library(ggeffects) #ggpredict
library(rstanarm) #stan models
library(shinystan) #stan model evaluation >launch_shinystan_demo()
library(loo) #loo() to compare fits between bayesian models
library(MuMIn) #dredge, model averaging
library(AICcmodavg) #aictab() for AICc model comparisons
library(RColorBrewer)
mc=tibble(read.csv("CB-mouthColor.csv", h=T))
nm=tibble(read.csv("nestlingMorph.csv", h=T))
nm[nm==0] <- NA
Morphometrics dataframe simplified for PCA
#subset id, billD, billW, TEC, skull, tarsus, and age
morphSimple.nm <- nm[,c(6,7,11,16:18,20,23)]
morphSimple.nm <- morphSimple.nm %>%
filter(between(age, 23, 27))
morphSimple.nm <- drop_na(morphSimple.nm)
#morphSimple.nm <- morphSimple.nm[-c(640,840,1186,1922,2196),]
morphForPCA <- as.matrix(morphSimple.nm[,c(3:7)])
Size PCA (same metrics as Townsend et al. 2010)
sizePCA <- prcomp(morphForPCA)
sizePCA.plot <- tibble(PC1=sizePCA$x[,1],
PC2=sizePCA$x[,2],
id=morphSimple.nm$id)
summary(sizePCA)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 4.0329 2.2578 1.7519 1.44284 0.58696
## Proportion of Variance 0.6056 0.1898 0.1143 0.07751 0.01283
## Cumulative Proportion 0.6056 0.7954 0.9097 0.98717 1.00000
sizePCA
## Standard deviations (1, .., p=5):
## [1] 4.0328741 2.2577553 1.7518941 1.4428431 0.5869643
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4 PC5
## tarsus -0.8008200 -0.08478806 0.15408891 -0.5716853 0.03051037
## billD -0.1357370 0.04076158 -0.06506740 0.1141950 -0.98114203
## billW -0.1200268 0.05105806 -0.98367877 -0.1007184 0.07223945
## TEC -0.4511211 0.68408742 0.04531459 0.5507967 0.15193321
## skull -0.3497801 -0.72150326 -0.04843046 0.5887395 0.09015093
weightByTarsusPlot <-morphSimple.nm %>%
ggplot(aes(x=tarsus, y=weight, label=id))+
geom_point(shape=1)+
geom_smooth()
ggplotly(weightByTarsusPlot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
massTarsus.mdl <- lm(weight~tarsus, data = morphSimple.nm)
plot(massTarsus.mdl)
morphSimple.nm$massTarsusResiduals <- residuals(massTarsus.mdl)
joinResiduals.df <- data.frame(id = morphSimple.nm$id,
residualsMassTarsus = morphSimple.nm$massTarsusResiduals)
new.mc <- left_join(mc, joinResiduals.df, by = "id")
mc$residualsMassTarsus
## Warning: Unknown or uninitialised column: `residualsMassTarsus`.
## NULL
Residuals (mass by tarsus) model for Saturation 1
#checked age * tarsus interaction...not significant
# sat1_resid.mdl <- lmer(sat1 ~ age + residualsMassTarus.x + sex + bodTemp1 + ambTemp1 + timeOut1 + (1|family),
# data = new.mc, REML = FALSE)
#
# summary(sat1_full.mdl)
Full model for Saturation 1
#checked age * tarsus interaction...not significant
sat1_full.mdl <- lmer(sat1 ~ weight + age + tarsus + sex + bodTemp1 + ambTemp1 + timeOut1 + (1|family),
data = mc, REML = FALSE)
summary(sat1_full.mdl)
## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: sat1 ~ weight + age + tarsus + sex + bodTemp1 + ambTemp1 + timeOut1 +
## (1 | family)
## Data: mc
##
## AIC BIC logLik deviance df.resid
## 1322.8 1356.7 -650.4 1300.8 151
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9805 -0.4147 0.1147 0.5441 1.7181
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 133.65 11.561
## Residual 94.32 9.712
## Number of obs: 162, groups: family, 88
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 190.47318 55.82130 3.412
## weight 0.09456 0.03766 2.511
## age -3.51493 0.62139 -5.657
## tarsus -1.34739 0.64571 -2.087
## sexF 5.74332 9.94909 0.577
## sexM 3.49846 9.74417 0.359
## bodTemp1 1.26644 1.42201 0.891
## ambTemp1 0.34515 0.40305 0.856
## timeOut1 0.01686 0.09879 0.171
##
## Correlation of Fixed Effects:
## (Intr) weight age tarsus sexF sexM bdTmp1 ambTm1
## weight 0.447
## age -0.157 -0.197
## tarsus -0.395 -0.748 0.075
## sexF -0.133 -0.006 0.006 0.121
## sexM -0.157 -0.029 0.004 0.082 0.973
## bodTemp1 -0.780 -0.114 -0.123 -0.162 -0.104 -0.046
## ambTemp1 0.261 -0.065 -0.016 0.209 -0.080 -0.093 -0.496
## timeOut1 -0.087 -0.028 -0.166 0.076 -0.013 -0.023 0.062 -0.065
Confidence intervals for Full Saturation 1 model
confint(sat1_full.mdl)
## Computing profile confidence intervals ...
## 2.5 % 97.5 %
## .sig01 8.89525628 14.58041653
## .sigma 8.29391209 11.52047360
## (Intercept) 80.19693221 301.08102026
## weight 0.01768514 0.17124532
## age -4.75215882 -2.28733611
## tarsus -2.62435838 -0.07139622
## sexF -14.06523964 25.62651294
## sexM -15.78560888 22.84723083
## bodTemp1 -1.58278201 4.08843893
## ambTemp1 -0.48547219 1.15254538
## timeOut1 -0.18884662 0.22891007
#plot model with sex
plot_model(sat1_full.mdl,
title = "Saturation 1")
#plot model without sex to view the smaller confidence intervals
plot_model(sat1_full.mdl,
title = "Saturation 1",
rm.terms = c("sex [F]","sex [M]"))
Reduced models: fixed effects removed based on confidence intervals from the full model
#sex random effect
sat1_rsex.mdl <- lmer(sat1 ~ weight + age + tarsus + bodTemp1 + ambTemp1 + (1|sex),
data = mc, REML = FALSE)
#remove timeOut
sat1_r1.mdl <- lmer(sat1 ~ weight + age + tarsus + sex + bodTemp1 + ambTemp1 + (1|family),
data = mc, REML = FALSE)
#remove ambTemp
sat1_r2.mdl <- lmer(sat1 ~ weight + age + tarsus + sex + bodTemp1 + (1|family),
data = mc, REML = FALSE)
#remove bodTemp
sat1_r3.mdl <- lmer(sat1 ~ weight + age + tarsus + sex + (1|family),
data = mc, REML = FALSE)
sat1_candset <- c(sat1_full.mdl,sat1_rsex.mdl,sat1_r1.mdl,sat1_r2.mdl,sat1_r3.mdl)
mdls <- c("full", "randSex", "r1", "r2", "r3")
aictab(sat1_candset, modnames = mdls, sort = TRUE)
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## full 11 1324.52 0.00 1 1 -650.38
## r1 10 1560.88 236.36 0 1 -769.84
## randSex 8 1608.09 283.57 0 1 -795.65
## r2 9 1810.46 485.94 0 1 -895.81
## r3 8 1831.02 506.49 0 1 -907.17
AICc model selection shows that the full model is better fit than any of the reduced models.
In summary, weight is positively associated with saturation 1. Age and tarsus are both negatively associated with saturation 1. Sex, ambient temp, and body temp are not significantly associated with saturation 1.